Strong electronic correlation in the Hydrogen chain: a variational Monte Carlo study 
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In this article, we report a fully ab initio variational Monte Carlo study of the linear, and periodic 
chain of Hydrogen atoms, a prototype system providing the simplest example of strong electronic 
correlation in low dimensions. In particular, we prove that numerical accuracy comparable to 
that of benchmark Density Matrix Renormalization Group calculations can be achieved by using a 
highly correlated Jastrow-antisymmetrized geminal power variational wave function. Furthermore, 
by using the so-called "modern theory of polarization" and by studying the spin-spin and dimer- 
dimer correlations functions, we have characterized in details the crossover between the weakly and 
strongly correlated regimes of this atomic chain. Our results show that variational Monte Carlo 
provides an accurate and flexible alternative to highly correlated methods of quantum chemistry 
which, at variance with these methods, can be also applied to a strongly correlated solid in low 
dimensions close to a crossover or a phase transition. 

PACS numbers: 71.30.-l-h, 73.21. Hb, 31.15.A-, 02.70.Ss 



I. INTRODUCTION 

The homogeneous (i.e., equispaced), linear, and peri- 
odic chain of Hydrogen atoms (hereafter, the H-chain) is 
commonly believed to be the simplest physical system de- 
scribed by the one-band, periodic, one-dimensional (ID) 
Hubbard Hamiltoniani (sec Eg. [TT|). This Hamiltonian 
is exactly solvable^ and its solution predicts the H-chain 
to be always a Mott- Hubbard (i.e., a strongly correlated) 
insulator. Indeed, it seems reasonable to model the H- 
chain by including only the Is orbitals and by neglecting 
the long-range tail of the Coulomb interaction, especially 
for large interatomic distances. As a consequence, the H- 
chain has been intensively studied to benchmark ab ini- 
tio approaches to strong electronic correlation^"— despite 
this atomic chain is not directly observable due to the 
well known Peierls' instability. 

Previous ab initio studies — mostly using methods of 
quantum chemistry — have been focused on finite (i.e., 
not periodic) linear chains, only. However, Periodic 
Boundary Conditions (PBC) analogous to the Born- 
von Karman boundary conditions used in solid state 
physics are better suited to investigate phase transitions 
or crossovers. Indeed, unwanted edge effects arc avoided 
by using PBC and a speed up of the convergence to 
the thermodynamic limit, i.e., the limit of infinite lin- 
ear chains, is expectediii^ Hence, an ab initio description 
of the low energy physics of a properly periodic H-chain 
is still missing. In this article we provide the first ex- 
haustive, fully ab initio^ variational description of peri- 
odic chains by using the same kind of variational wave 
function for both the weakly and the strongly correlated 



regimes, i.e., for both small and large interatomic dis- 
tances. 

From previous studies, it is known that in the strongly 
correlated regime — i.e., for interatomic distances, a, 
larger than a certain critical distance, Oc — the ground 
state of the finite H-chains is characterized by a huge 
degeneracy of the natural orbitals^ which leads to an 
almost uniform natural orbital population, narrowly dis- 
persed around This behavior has to be contrasted 
with the weakly correlated regime (a < Oc) for which al- 
ready the restricted Hartree-Fock reference — that yields 
either doubly occupied of empty natural orbital — is quite 
accurateJ^ 

In many cases part of the static correlation which 
characterizes the strongly correlated regime can be ef- 
fectively recovered by means of an unrestricted Hartree- 
Fock calculation, or cquivalently, by means of a spin- 
polarized density functional theory calculation within the 
local density approximationji Although justified for finite 
systems, a spin-polarized approach implies a mean-field 
antiferromagnetic order, which cannot be trusted in the 
thermodynamic limit, because true antiferromagnetism 
is not possible for ID solids.— 

The Density Matrix Renormalization Group (DMRG) 
provides an, in principle exact, algorithm to compute the 
electronic structure of ID and almost- ID systems, al- 
though in practice limited by the size of the orbital basis 
seti^ It works very efficiently also when other highly cor- 
related methods fail, e.g., configuration interactioni^ is 
not applicable if the system size is too large and the stan- 
dard coupled cluster singles and doubles plus perturba- 
tivc triples becomes unstable in ID for large interatomic 
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distances 4 Nevertheless, even the very favorable numeri- 
cal efficiency of DMRG is lost for a gapless (i.e., metallic) 
chain. In this case, also the MoUer-Plesset perturbative 
approach is not straightforwardly applicable due to the 
numerical issues triggered by the vanishing small gap. 

II. COMPUTATIONAL METHODS 

Among the possible alternatives to standard quan- 
tum chemical approaches?^ one can consider non- 
perturbative quantum Monte Carlo (QMC) methods like 
VMC or the, in principle more accurate. Diffusion Monte 
Carlo (DMC) In fact, a direct application of DMC to 
a homogeneous chain raises some ergodicity issues when 
the electrons are very localized about the nuclei, i.e., in 
the strongly correlated regime As a consequence, to 
date only alternating (i.e., dimerized) chains have been 
investigated by means of DMGi^ On the other hand, 
VMC can be made ergodic by tailoring the sampling and 
by improving the variational many-body wave function, 
so that it remains very effective also close to crossovers 
and phase transitionsjii^ 

There have been dramatic theoretical advances in the 
quality of the variational wave function during the last 
decade, so that it is now possible to achieve chemical 
accuracy for atoms, ions, small molecules and even peri- 
odic systemsiii"— These quantitative improvements have 
been made possible by new stochastic optimization tech- 
niques which can optimize variational wave functions de- 
pending on hundreds free parameters 

In our variational calculations, we have used a Jastrow- 
antisymmctrized geminal power ( JAGP) variational wave 
function for (even) interacting electrons fi^ii^ 

N/2 

^t,{R)^J{R)AX{W],^.), (1) 

i=l 

where R ~ |r*[, . . . , rj^^j' ''i ' ■ • • i ^72} ^A^- 
dimensional coordinate vector, A is the antisymmetriza- 
tion operator, ^{x^ ,y^) is a symmetric function describ- 
ing a singlet electron pair and 

N 

J{R) = exp^[7.(r,,,) + fin, r,)] , (2) 

is the Jastrow factor. 

By using Eq. [1] one can accurately describe both static 
and dynamics correlations. Indeed, the antisymmetrized 
geminal power, 

'4n»=^i^*(^I,^) provides a very com- 
pact multideterminant reference, while the Jastrow fac- 
tor takes into account the dynamic correlation from; i) 
A short-range homogeneous electron-electron interaction 
through the term u{rij) which just depends on the dis- 
tance, j = 1?^; — rj\, between paired electronsj^ ii) The 
inhomogeneous electron-electron-nucleus and electron- 
electron-nucleus-nucleus interactions through the term 



f{fi,fj), which depends separately on the coordinates of 
the paired electrons, and it can also describe long-range 
electronic correlations. 

Results showed in this article have been obtained us- 
ing a short-range homogeneous Jastrow factor, u{rij) = 
(&/2)(l — e~'""'-j). In addition, to asses the sensitivity 
of these results on the detail of the wave function, we 
have also considered a long-range homogeneous Jastrow 
factor, w(r,;j) = rij/2{l + brij) (see Fig. S]). In both 
cases, 6 is a variational parameter, and the electronic 
cusp conditions'^ are automatically satisfied. 

The functions ^{x^,y^) in Eq. [T]and f{x,y) in Eq. [5] 
are expanded using (in principle different) nonorthogonal 
atomic orbitals?^ {</),} and {ifi} so that 

$(f^,/) = ^ A„,;3<Pa(x^)</'^i(y^) , 

a, 13 

f{x,y) =^ga,l3<l)a{x)(l)/3{y) ■ (3) 

a,l3 

In particular, up to 3s orbitals for the geminal part and 
2s2p for the Jastrow part have been considered in this 
workj^ 

In principle, all the entries of the matrices Aq.^ and 
ga.p in Eq. [3]are variational parameters to be optimized. 
However, by taking advantage of the symmetries of the 
periodic linear chain?^ and by using an alternative, mini- 
mal expansion in terms of molecular orbitals the actual 
number of independent variational parameters to opti- 
mize is reduced, so that the optimization can be effec- 
tively performed by the method describes in Ref. [l^. 

All the VMC calculations have been performed by us- 
ing the TurboRVB package,— starting from a density 
functional theory (in local density approximation) calcu- 
lation employing the orbital basis set, of the gem- 
inal part (see Eq. [S]). This preliminary step is done to 
speed up the convergence of the following VMC optimiza- 
tion, while avoiding an uncontrolled bias. 

Finally, we have used a supercell with PBGi^ in all 
three Cartesian directions to model the periodic linear 
chain. To avoid unphysical self-interaction of the chain 
with its periodic replicas, the transverse dimensions of 
the supercell have been taken as min(16a, 80) a.u., where 
a is the interatomic distance. According to the supercell 
formalism, the thermodynamic limit, 00, can be 

achieved by increasing the number. . of H atoms in the 
supercell. In particular, where not otherwise indicated, 
by Hjv we mean a periodic chain with A'^ atoms in the 
supercell. 

To identify the weakly and the strongly correlated 
regimes, we have used the so-called "modern theory of 
polarization" .— This theory also provides a way to dis- 
criminate between a metal and an insulator alternative 
to the knowledge of the (many-body) charge gap, which 
in fact is not accessible by a variational ground state cal- 
culation. 

In practice, by VMC one computes the expectation 
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values of the complex polarization, zjq^ 

zjv = (*iv|e**^»'-"|*Ar), (4) 

where r| is the component of parallel to the chain axis. 
Then the electronic localization length, Aat, is obtained 
as 



In \zm\ 



N 



(5) 



where L is the longitudinal dimension of the supercell 
and N the number of H atoms in the supercell. 

From previous studies,— i^i^ one expects a huge degen- 
eracy of the natural orbitals when the electrons are very 
localized about the nuclei, i.e., when Xn/cl <C 1. Be- 
sides, the theory says that, in the thermodynamic limit, 
iV — >^ 00, a metal is characterized by a vanishing modu- 
lus of the complex polarization, |zjv| — J* 0, while in the 
insulating case |zjv| — > l"^*^ 



III. RESULTS 
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FIG. 1. (Color online) Panel (a): Total energy per atom as a 
function of the interatomic distance from VMC calculations 
of periodic chains with 18, 34, 50, and 66 H atoms in the 
supercell. [Data are almost superimposed at the scale of this 
figure, see also Tab. ID] Inset (b): Comparison between the 
total energy per atom of a finite Hsu chain obtained by VMC 
and DMRG^. 

In Fig. [TJa), we show the convergence of the total en- 
ergy per atom by increasing the number of H atoms per 
supercell for several interatomic distances. We note that 
the H50 periodic H-chain is already well converged at the 
scale of this figure. To follow the fine detail of the con- 
vergence, the values of the total energy per atom details 
have been also listed in Tab. HI 

In Fig. [TJb), a direct comparison between the VMC 
total energy for the H50 finite chain and the benchmark 
DMRG results obtained by using a ST0-6G basis seli^ 
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TABLE I. Total energy per atom as a function of the inter- 
atomic distance, a, for the same periodic chains of Fig. [TJa). 
The VMC error on the last digit is indicated in parenthesis. 



demonstrates the accuracy of our optimized JAGP vari- 
ational wave functioui^ In this case, to have a fair com- 
parison against the DMRG data, PBC have not been em- 
ployed to obtain the VMC results showed in Fig. [TJb). 
The difference between the total energy of H50 chains 
with and without PBC and the same interatomic dis- 
tance is of the order of few mHa per atom. 

Having verified the quality of the variational wave 
function, in Fig. [IJa) we plot the electronic localization 
length, Atv, in units of the interatomic distance, a, as a 
function of a. For all the supercells considered, we find 
that 



- if a < flc 
i^^ if a > Ur 



(6) 



where 77 ~ 0.5 and Oc — 3.5 a.u. This critical behavior is 
also in agreement with the sudden switch from \z\ ~ to 
\z\ ~ 1 visible in Fig. [5l^b), i.e., to the crossover between 
a (finite-size) metal and an insulator, namely a Mott- 
Hubbard insulator.— 

To further characterize the nature of the weakly and 
strongly correlated regimes of the H-chain, we have in- 
vestigated the spin-spin, 

fU^-j) = {^N\Si'^Si^^\^N) (7) 
and the dimer-dimer 

Ud{^-J) = {^N\s^'^si^+'^si^^si^+'M^N) (8) 

correlation functions, where S'i*^ measures the transverse 
component of the electronic spin about the ith H atom 
of the chain. By neglecting logarithmic corrections, we 
have fitted these functions by22, 
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FIG. 2. (Color online) Panel (a): Electronic localization 
length, Ajv, divided by the interatomic distance, a, as a func- 
tion of a, for the same chains of Fig. [TJa). Panel (d): modulus 
of the complex polarization, |2jv| as a function of the inter- 
atomic distance for the same chains of Panel (c). 

for the spin-spin and the dimer-dimer cases, respectively. 
The parameters, a^^, bss, Kss, add, bdd, and Kdd have 
been fitted independently for each value of the inter- 
atomic distance, a, and the number of H atoms in the 
supercell. 

Results for the scaling exponents, Kss and Kss, of the 
H-chain arc reported in Fig. ^a) and Fig. ^c) , respec- 
tively. For comparison, in Fig.[3l^b) and Fig.|3l^d) we show 
the scaling exponents obtained by solving numerically^^ 
the iV-site Hubbard model Hamiltonia»i with PBC (i.e., 
the simulation cell is folded so that the {j+N)th and jth 
sites represent the same atom) 

N N 
j=l a=t.i j = l 

(11) 

with a number of sites, N, correspondent to the 
number of H atoms in the chain supercell. In 
Eq. [TTl the creation(annihilation) operator c] g.(c^ g.) cre- 
ates (annihilates) an electron of spin a at site j, while 
rij^a- = c]^c^ g.. Since one expects (for U fixed) that 
\og{U/t) oc a, we have shown the Hubbard exponents 
as a function of U jt using a semi-log plot. 

The behavior of the scaling exponent Kss is very sim- 
ilar in the two cases, i.e., H-chain in Fig. EJa) and Hub- 



FIG. 3. (Color online) Panel (a): Scaling exponent Kss of 
the spin-spin correlation function (see Eq. [Qj as a function 
of the interatomic distance, for same H-chains of Fig. [T]and 
Fig. [2] Panel (b): Scaling exponent K^a of the Hubbard model 
with a corresponding number of sites (see text). [Data are 
superimposed at the scale of this figure.] Panel (c) and (d): 
same as panel (a) and (b), but for the dimer-dimer correlation 
function (see Eq. llOp. 

bard model in Fig. El^b). However, small but notice- 
able discrepancies between the H-chain and the Hubbard 
model for the scaling exponent Kdd arc found by com- 
paring Fig. |3l[c) and Fig. [3ld). 

For large interatomic distance a, both the H-chain ex- 
ponents behave as expected for the Hubbard model (ne- 
glecting logarithmic corrections) i.e., Kdd ^ Kss ^ 1^ 
Less conclusive results can be inferred from the weakly 
correlated regime for a < Oc- Finite-size effects are re- 
sponsible for the deviation of the scaling exponents from 
their thermodynamic values in the case of the Hubbard 
model. The same effects also mask possible discrepancies 
in the thermodynamic limit of the H-chain and the Hub- 
bard model. Further numerical investigations are needed 
to provide conclusive results on a possible metal-insulator 
transition at a finite value, Uc, of the interatomic distance 
of the H-chain. (see Sec. ITV)) . 

Finally, we investigate in more detail the capability 
of the variational JAGP wave function to describe the 
Mott- Hubbard insulating phase of the H-chain for a > Qc- 
In particular, wc focus on the a = 5.0 a.u. case and 
we consider some variants of the JAGP variational wave 
function, Eq. [T] 
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FIG. 4. (Color online) Panel (a): Localization length as a 
function of the inverse of the number of H atoms in the su- 
percell for different parametrizations of the JAGP variational 
wave function (see text). [Data 'sr-hJ+iJ' and 'Ir-hJ+iJ' are 
superimposed at the scale of this figure, see also Tab.HIl] Inset 
(b), the corresponding total energy per atom. 



In Fig. m^a) we plot the electronic localization length, 
Xn, obtained by optimizing different JAGP variational 
wave functions. We consider the following cases (for the 
notation, see previous Section): i) sr-hJ+iJ, correspond- 
ing to Eq. [1] with short-range homogeneous Jastrow fac- 
tor plus the inhomogeneous part. This is the standard 
case considered elsewhere in this article; ii) lr-hJ-|-iJ, as 
the sr-hJ-t-iJ wave function, but with a long-range homo- 
geneous Jastrow factor, instead; iii) sr-hJ, and iv) Ir-hJ, 
which are as the sr-hJ+iJ and Ir-hJ-fiJ wave functions, 
respectively, but without the inhomogeneous Jastrow fac- 
tor, i.e., with f{ri,rj) = 0. 

We find that the localization length, Aat, is well con- 
verged at the scale of Fig.SJJa) if the inhomogeneous Jas- 
trow factor is included, regardless of the choice of the 
homogeneous part. If this factor is not included, values 
of Ajv almost twice as large are found and they slightly 
increase with the system size, showing that the homoge- 
neous Jastrow alone is not enough to give an accurate 
description of the Mott-Hubbard insulating phase. Val- 
ues of the total energy per atom reported in Fig.l^Jb) also 
confirm the relevance of the long-range, inhomogeneous 
Jastrow factor to improve the accuracy of the VMC de- 
scription of a Mott-Hubbard insulator. To follow the fine 
detail of the convergence, the values of the total energy 
per atom have been also listed in Tab. [TTl 

Our findings are in agreement with previous varia- 
tional studies of lattice models with short-range inter- 
actions which showed that a correct description of the 
Mott-Hubbard insulating phase can be achieved only by 
combining a Gutzwiller projector and a long-range Jas- 
trow factor i^i^ 

In fact, although the single H atoms are on average 
neutral, charge fluctuations which give, e.g., virtual H''"- 
H~ pairs, are always possible. Such charge fluctuations 
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TABLE II. Total energy per atom as a function of the num- 
ber, A'^, of H atoms in the supercell, for the same chains of 
Fig. m^b). The VMC error on the last digit is indicated in 
parenthesis. 



are in fact strongly suppressed at large interatomic dis- 
tance, a (and, in the Hubbard model, for large U). There- 
fore, at large a, in order to prevent an instability of the 
Mott-Hubbard insulator toward a metallic state driven 
by the quantum fluctuations f2i the H'''-H~ pairs must 



be bound. In the context of the Hubbard model 
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has been demonstrated that such binding mechanism can 
be included in the variational wave function by means of 
an inhomogeneous Jastrow factor, analogous to the the 
second term in Eq. [2] Crucially, the matrix element ga^p 
in the expansion of the inhomogeneous Jastrow factor 
(See. Eq.|31) can be nonvanishing also for pairs of orbitals 
(labeled by the Greek indices) which are far apart. In- 
deed, such long range correlation is necessary to bind the 
virtual H+-H~ pairs and, along with the homogeneous 
Jastrow factor, provides an accurate way to model elec- 
tron localization in the Mott-Hubbard insulating phase, 
as shown in Fig. IH^a). 



IV. DISCUSSION AND CONCLUSIONS 

In this article, we have investigated the ground state 
of a homogeneous, linear, and periodic chain of Hydro- 
gen atoms (or H-chain) from first principles, by means of 
a state-of-the-art variational Monte Carlo approach. In 
fact, using a highly correlated Jastrow-antisymmetrized 
geminal power variational wave function allowed us to 
obtain a total energy per atom comparable to bench- 
mark density matrix renormalization group calculations 
Furthermore, by using the so-called "modern theory of 
polarization" ^ we have characterized the crossover be- 
tween the weakly correlated (for small interatomic dis- 
tances) and the strongly correlated regimes (for large in- 
teratomic distances) through a single, simple parameter, 
i.e., the electronic localization length. Our results ex- 
tend to the properly periodic H-chain previous results ob- 
tained by studying long, yet finite, chainsj^i^^ In partic- 
ular, we confirmed that the crossover between the weakly 
and strongly correlated regimes of the H-chain corre- 
sponds physically to a crossover between a (finite-size) 
metallic and an insulating phase. Finally, by studying 
the asymptotic behavior of the spin-spin and dimcr-dimcr 
correlation functions, we have also verified that the insu- 
lating phase of the real H-chain is of the Mott-Hubbard 
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type, as expected from the Hubbard modelji Since the 
correct description of such a correlated insulator is be- 
yond the possibility of density functional theory in any of 
its conventional local or semi-local approximations, one 
can think of using our findings to devise an improved 
class of functionals. We are currently exploring this pos- 
sibility. 

Intriguingly, we have found small but noticeable devia- 
tions from the behavior predicted by the Hubbard model 
in the case of the scaling exponent of the dimer-dimer 
correlation function predicted by the Hubbard model.— 
These deviations can be possibly due to the finite-size 
scaling or to a true discrepancy between the H-chain and 
the ID Hubbard model in the thermodynamic limit. 

It will be interesting to check for possible new low en- 
ergy physics of the H-chain at variance of the one-band, 
ID Hubbard model predictions. These might be origi- 
nated by: i) The long-range Coulomb repulsion — indeed 
inefficiently screened in ID systems — not included in the 
Hubbard model; ii) The atomic orbital polarization — 
essential to describe non-covalent contribution to the 
bonding — not representable in terms of Is orbitals, only. 
In particular, relative simple elaborations of the Hubbard 
model, which just contain next-nearest neighbor interac- 
tion, already predict a rich phase diagram even for a ID 
system.— Besides, it is known that long-range Coulomb 
repulsion can yield gapless charge excitations (plasmons) 
in ID, as observed in experimentsi^ 

Of course, more accurate finite-size extrapolation is 
desirable, although not possible with our current com- 
putational resources. In particular, the use of diffusion 
Monte Carlo to improve the variational results might be 



also beneficial, but in practice still highly problematic 
due to the well known ergodicity issues in dealing with 
strong electronic localization in ID systems.— 

In conclusion, given that the homogeneous, linear and 
periodic chain of Hydrogen atoms is becoming a stan- 
dard test-case for ab initio approaches to strong elec- 
tronic correlation,i^i^ the results reported in this arti- 
cle show that variational Monte Carlo (with a highly 
correlated Jastrow-antisymmetrized geminal power vari- 
ational wave function) can provide an accurate and flex- 
ible alternative to highly correlated methods of quantum 
chemistry. Besides, and at variance with most of the 
methods of quantum chemistry, variational Monte Carlo 
can be successfully employed to study a strongly corre- 
lated solid in low dimensions close to a crossover or a 
phase transition!^ 
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